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Strong coupling of electronic and vibrational degrees of freedom entails a low-bias suppression of 
the current through single-molecule devices, termed Franck-Condon blockade. In the limit of slow 
vibrational relaxation, transport in the Franck- Condon-blockade regime proceeds via avalanches of 
large numbers of electrons, which are interrupted by long waiting times without electron transfer. 
The avalanches consist of smaller avalanches, leading to a self-similar hierarchy which terminates 
once the number of transferred electrons per avalanche becomes of the order of unity. Experimental 
signatures of self-similar avalanche transport are strongly enhanced current (shot) noise, as expressed 
by giant Fano factors, and a power-law noise spectrum. We develop a theory of the Franck- Condon- 
blockade regime with particular emphasis on effects of electron cotunneling through highly excited 
vibrational states. As opposed to the exponential suppression of sequential tunneling rates for 
low-lying vibrational states, cotunneling rates suffer only a power-law suppression. This leads to a 
regime where cotunneling dominates the current for any gate voltage. Including cotunneling within a 
rate-equation approach to transport, we find that both the Franck-Condon blockade and self-similar 
avalanche transport remain intact in this regime. We predict that cotunneling leads to absorption- 
induced vibrational sidebands in the Coulomb-blockaded regime as well as intrinsic telegraph noise 
near the charge degeneracy point. 

PACS numbers: 73.63.-b, 81.07.Nb, 72.70.+m, 63.22.-fm 
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I. INTRODUCTION 

According to many recent experiments,^'^ transport 
through single-molecule junctions is characterized by 
Coulomb-blockade behavior and exhibits a Kondo effect 
at low temperatures, in close analogy with experiments 
on quantum dots. An important distinction between 
transport through quantum dots and single molecules lies 
in the coupling to phonon degrees of freedom. While 
charge carriers typically interact with a continuum of 
phonon modes in quantum dots, molecules are character- 
ized by a discrete spectrum of vibrational modes. Recent 
experiments on single-molecule junctions have shown 
that these discrete excitations are reflected in current- 
voltage characteristics {IVs) as vibrational sidebands. "^"^ 
This experimental effort has motivated much theoretical 
work.'^-i" It is now understood that transport through 
single-molecule junctions exhibits a number of qualita- 
tively different regimes, depending on parameters such 
as the strength of the electron-phonon coupling and the 
vibrational relaxation rate, see e.g. Refs. 11-14. 

The importance of Coulomb-blockade physics in single- 
molecule junctions suggests that molecule and leads are 
often weakly coupled from the point of view of electronic 
transport. Hence, perturbation theory in the tunneling 
amplitude between molecule and leads should be good 
starting point for a theoretical description. Such a per- 
turbative treatment is valid if the tunneling-induced level 
width r is small compared to the temperature T. Then, 
the lowest order in this expansion results in sequential 
tunneling, which corresponds to the transfer of one elec- 
tron from a lead onto the molecule or vice versa. Sequen- 
tial tunneling is dominant as long as the electronic level 



is situated between the Fermi energies of the leads. By 
contrast, if the level is located outside of the bias win- 
dow, sequential tunneling is exponentially suppressed, 
and (next-to-leading order) cotunneling processes take 
over.^^ This marks the onset of Coulomb blockade. In 
many experiments, a gate electrode can be used to tune 
the system between the sequential tunneling and the co- 
tunneling regime. 

The coupling of electron tunneling to molecular vibra- 
tions adds a wealth of new physics. This coupling origi- 
nates from the fact that the nuclear configuration of the 
molecule needs to adjust to the addition of electrons to 
the molecule. Thus, electron transfer is generally accom- 
panied by excitation (or deexcitation) of molecular vibra- 
tions. In the simplest case, the electron-phonon coupling 
can be modeled by a charge-dependent displacement of 
the molecular potential surface. 

In two recent papers, ^^'^^ two of us have investigated 
the regime of strong electron-phonon coupling, where the 
displacements of the potential surfaces are large com- 
pared to the quantum fluctuations of the nuclear config- 
uration in the vibrational ground state. As a result, the 
overlap between low-lying vibrational states is drastically 
suppressed. It is this suppression of the Franck-Condon 
matrix elements which leads to a strong low-bias suppres- 
sion of the sequential-tunneling current, which we termed 
Franck-Condon (FC) blockade. It is a hallmark of the FC 
blockade that it cannot be lifted by means of a gate volt- 
age, in striking contrast with other blockades such as the 
Coulomb blockade. 

It is evident that the FC blockade is more severe for 
molecules with fast vibrational relaxation (arising, e.g.. 
due to interactions with the substrate, radiation, or cou- 



pling to other vibrational modes) . The reason is that for 
slow vibrational relaxation, the transport current excites 
the molecular vibrations, and transport involves excited 
vibrational states for which the suppression of the FC 
matrix elements is less pronounced. Nevertheless, the 
Franck-Condon blockade in the current-voltage charac- 
teristics exists irrespective of the strength of vibrational 
relaxation. 

It is thus remarkable that the transport mechanisms in 
the limits of weak and strong relaxation are in fact fun- 
damentally different, as revealed by current shot noise. 
In the sequential-tunneling approximation, strong relax- 
ation leads to purely sub-Poissonian current noise, as 
commonly expected for transport in fcrmionic systems. 
Electrons are transferred one by one, slightly avoid- 
ing each other due to the Pauli principle. By con- 
trast, for weak vibrational relaxation transport is domi- 
nated by self-similar avalanches of electrons with giant 
Fano factors reflecting the large number of electrons per 
avalanche. The occurence of avalanches is explained by 
the exponential growth of sequential-tunneling rates with 
increasing vibrational quantum number. Whenever the 
molecule is in its vibrational ground state, the suppressed 
FC matrix elements require a long waiting time before an 
electron tunnels. Most likely, the transition involves an 
excitation of the vibrations. Now, having left the vibra- 
tional ground state, the accessible transitions have faster 
rates, leading to rapid tunneling dynamics and avalanche 
formation. If, eventually, the system undergoes a transi- 
tion back into the vibrational ground state, the avalanche 
terminates and another long waiting time is induced. 
The self-similarity arises, since all previous arguments 
can also be applied to the first excited vibrational state. 
As an experimental fingerprint of the self-similarity, we 
have predicted a characteristic power-law decay of the 
current noise as a function of frequency. 

In this paper, we investigate the cotunneling contri- 
butions to current and noise in the FC blockade regime. 
The inclusion of such higher-order processes is motivated 
by the insight that sequential-tunneling rates are expo- 
nentially suppressed in the electron-phonon coupling. By 
contrast, cotunneling may benefit from the participation 
of highly excited virtual phonon states with much better 
overlap, see Fig. 1. In this way, the relevant FC matrix 
elements are no longer exponentially suppressed, which 
may render cotunneling relevant despite its suppression 
in powers of the small quantity T/ksT. The central ques- 
tion which we address in this paper is whether and if 
so, in which way, cotunneling affects the phenomenology 
of the Franck-Condon blocakade as sketched above. We 
investigate this question within the rate-equations ap- 
proach, extending the formalism to include cotunneling 
corrections. 

For the convenience of readers who are less interested 
in the technical details, we begin our presentation in Sec- 
tion II with a summary of our principal results. In Sec- 
tion III we present the model that our analysis is based 
upon. The computational formalism and the question 
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FIG. 1: (Color online) (a) Sequential tunneling versus (b) 
cotunneling for strong electron-phonon coupling. In the FC 
blockade regime, sequential processes as depicted in (a) are 
strongly suppressed due to the exponentially small overlap of 
harmonic oscillator wavefunctions. For cotunneling processes 
(b), this suppression is partially lifted due to contributions 
from highly excited virtual phonon states. 



of regularized cotunneling rates are addressed in Section 
IV. In particular, we discuss the cotunneling rates in the 
FC blockade regime, and determine their scaling with the 
electron-phonon coupling. Section V contains the inves- 
tigation of cotunneling corrections to the linear conduc- 
tance and nonlinear current-voltage characteristics. Fur- 
ther insight into the transport mechanism and the fate 
of avalanches is provided in Section VI, where we discuss 
the current shot noise beyond the sequential-tunneling 
approximation. Finally, the conclusions from this anal- 
ysis are summarized in Section VII. Technical details 
concerning the calculations are deferred to appendices. 



II. QUALITATIVE DISCUSSION 

A. The fate of the FC blockade 

The increased overlap of vibrational wavefunctions ex- 
ploited by cotunneling processes indeed makes cotunnel- 
ing corrections significant in the FC blockade regime. In 
contrast to sequential-tunneling rates, we find that co- 
tunneling rates are not exponentially, but rather alge- 
braically suppressed in the electron-phonon coupling, see 
Eq. (14). Consequently, for sufficiently strong electron- 
phonon coupling, cotunneling dominates the linear con- 
ductance over the entire gate-voltage range, see Fig. 4. 
In particular, the paradigm of a gate-controlled crossover 
between sequential tunneling and cotunneling breaks 
down in this regime, and the linear conductance devel- 
ops a kink at the degeneracy point. Cotunneling also 
significantly increases the low-bias current well beyond 
the linear-response regime, see Fig. 6. More importantly 
however, we find that the remaining strong algebraic sup- 
pression of cotunneling rates still leads to a current reduc- 
tion. Therefore, a central finding is that the FC blockade 
as such remains intact. 
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B. Persistence of avalanches and enhanced Fano 
factors 

Another central result is that the avalanche dynam- 
ics, found in the sequential-tunneling approximation for 
weak vibrational relaxation, persists in the presence of 
cotunneling. This can be understood by noting that co- 
tunneling is dominant, whenever the molecule is close 
to its vibrational ground state. However, the avalanches 
involve excited vibrational states for which cotunneling 
is in fact a small perturbation. Nevertheless, the long 
waiting times between avalanches are now governed by 
cotunneling, leading to Fano factors which are smaller 
than predicted by the sequential-tunneling approxima- 
tion, although still strongly enhanced compared to the 
Poissonian limit, see Fig. 9. In a similar vein, the power- 
law behavior of the noise power spectrum remains valid, 
but is reduced in the frequency range over which it ap- 
plies, see Fig. 10. 

C. Absorption-induced vibrational sidebands 

For weak vibrational relaxation, cotunneling leads to a 
striking new feature: the appearance of additional vi- 
brational sidebands inside the Coulomb-blockaded re- 
gion, see Fig. 7. Inside the Coulomb-blockaded regime, 
molecular vibrations can be excited by inelastic cotun- 
neling. Once molecular vibrations arc excited, sequential 
tunneling processes become possible when accompanied 
by absoption of a phonon. This leads to new vibra- 
tional sidebands which were in fact observed in exper- 
iments on suspended carbon nanotubes.^® Interestingly, 
such absorption-induced sidebands are expected to be a 
generic feature which is not limited to strong electron- 
phonon coupling. 

D. Inherent telegraph noise 

For strong vibrational relaxation, avalanche dynamics 
is ruled out and current noise remains sub-Poissonian 
in the sequential-tunneling approximation. Surprisingly, 
we find that cotunneling causes a strong enhancement of 
current noise at low bias due to inherent telegraph noise, 
see Figs. 2 and 8. Traces of this telegraph behavior are 
also found for weak vibrational relaxation, rendering this 
effect robust with respect to the relaxation rate of the 
system. Its origin is the unusual situation that cotun- 
neling rates are larger than sequential-tunneling rates, 
Wcot ^ W^seq- Thus, Sequential tunneling leads to a slow 
switching between the two charge states n = and 1, 
which modulates the fast cotunneling dynamics. The co- 
tunneling rates in these two states are different in gen- 
eral, which leads to two distinct current states. This is 
reflected in large Fano factors F ~ AWcot/W^seq ^ 1; 
where AWcot denotes the difference of cotunneling rates 
in the two charge states. 
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FIG. 2: Coarse-grained current I{t) as obtained by Monte- 
Carlo simulation for equilibrated phonons at strong electron- 
phonon coupling (A = 5), with bias and gate voltage fixed 
to eV = O.Shu) and = 0.5?lu;, respectively. The time- 
dependence of the current clearly reveals the current telegraph 
noise caused by a slow switching between the two molecu- 
lar charge states via sequential tunneling. (fcsT — 0.05/ta;, 
Ta = 0.02?ia;; coarse graining over 300 events) 

III. MODEL 

The Anderson-Holstein Hamiltonian constitutes the 

minimal model which captures Coulomb blockade, the 
Kondo effect, as well as vibrational sidebands. 
Several assumptions and simplifications enter into this 
model: (i) Transport is assumed to be caused by tunnel- 
ing through one spin-degenerate orbital of the molecule, 
(ii) One mode of molecular vibrations is taken into ac- 
count within the harmonic approximation. These ingre- 
dients result in the Hamiltonian H = ifmoi + ^^leads + Ht 
with 

-ffmoi =ednd + UndfUdi + hojb^b + Xfuj{b^ + b)nd, (1) 

describing the electronic and vibrational degrees of free- 
dom of the molecule, 

-ffleads — ^ ^ ^ ^ (^ap ^ Ma ) f-apcr ^apcr (2) 
a=L,R p,<7 

the noninteracting leads, and 

Ht^ J2 E (^-pcipc.'^- + h-c.) (3) 

a=L,R p,<7 

the tunneling between leads and molecule, (iii) In addi- 
tion, electronic relaxation is taken to provide the shortest 
time scale in the problem, so that the metallic electrodes 
are described by Fermi seas in thermal equilibrium at all 
times. 

The validity and various extensions of this model in 

the context of transport through single molecules have 
been discussed, e.g., in Refs. 9,10,12,24-26. However, 
the inclusion of generalizations would only obscure the 
central points of the present paper, and we therefore re- 
strict ourselves to the model constituted by Eqs. (l)-(3). 



4 



For simplicity, we consider the case of symmetric volt- 
age splitting, i.e. the bias voltage V enters the left and 
right Fermi energies as ^iL,R = ±eV/2. (Note that all en- 
ergies are measured with respect to the zero-bias Fermi 
energy.) All results are derived for the wide-band limit 
with constant densities of states pl,r and momentum- 
independent tunneling matrix elements ta- 

There are four different energy scales relevant to the 
analysis of transport regimes of the Anclerson-Holstoin 
model: the phonon energy hv, the tunneling- induced 
level width F = '^T'^a Pa\taf , and the energy scales 
set by the Kondo temperature Tk and the ambient 
temperature T. Here, we are interested in the "high- 
temperature" limit of the model, i.e. we study the regime 
where ksTK, F ^ huj, ksT. In this case, the tunneling 
Hamiltonian Hj^ may be treated perturbatively, and it 
has been shown that the equation of motion for the re- 
duced density matrix of the molecule reduces to simple 
rate equations. 

The electron-phonon coupling, parametrized by the di- 
mensionless coupling parameter A, can be treated exactly 
by means of the Lang-Firsov canonical transformation.^^ 
This transformation amounts to rcnormalizations of the 
single-particle energy —^ e'^ = Cd — X^hui, and of the 
charging energy U ^ U' = U — 2}?tkJ (polaron shift), 
and it introduces a translation operator into the tun- 
neling matrix elements, ta exp[— A(6t — b)]ta- This 
shift ^ A in the normal coordinate of the harmonic os- 
cillator precisely corresponds to the fact that the equi- 
librium distances between the nuclei of the molecule 
change whenever an additional electron is present in the 
molecule. The transformed Hamiltonian reads H' = 

^mol + -f^leads + ^^T ^ith H(^^^^ = illeads, and 

-f^moi + U'ud^Udi + hwb^b, (4) 

= E E {ta^e-''-'''-'^cl^„d„ + h.c.) . (5) 

a=L,rt, p.CT 

In the following, we will proceed with the transformed 
Hamiltonian and drop all primes. 

The system is expected to display qualitatively dif- 
ferent behavior depending on the vibrational relaxation 
rate as well as the strength of the electron-phonon cou- 
pling. Throughout this paper, we will focus on the case 
of strong electron-phonon coupling, A 3> 1,'^^ and non- 
negative effective charging energies, U — 2\^hjj > 0. The 
latter inequality is always satisfied for sufiiciently strong 
Coulomb repulsion between electrons inside the molec- 
ular orbital. Unless explicitly stated otherwise, we will 
concentrate on the limit U ^ oo, applicable whenever U 
is large compared to all other relevant energy scales. 



IV. FORMALISM 

In this section, we discuss the formalism which consti- 
tutes the basis for our analysis. The main input for the 



transport calculations consists of the transition rates in- 
duced by electron tunneling. In the following subsection, 
we discuss how to obtain these rates, especially in the 
case of cotunneling where the naive application of Fermi's 
golden rule leads to a divergence. The subsequent subsec- 
tion explains the rate-equations approach for calculating 
the steady-state current as well as the noise spectrum. 



A. Transition rates 

Regarding the tunneling Hamiltonian Ht as a pertur- 
bation, one can proceed to calculate rates for transitions 
N ) ^ I / ) by an expansion of the T-matrix and applying 
Fermi's golden rule, 

j,f = ^\{f\T\t}\'6iEi-Ef) (6) 
= ^\{f\HT + HtGoHt + ■■■ \ i)f6{Ei-Ef). 

Here, Gq = [Ei - H^^i - -f^ieads + is the free re- 

tarded Green's function. Generally, the initial and fi- 
nal states I i ), | / ) still involve degrees of freedom of the 
leads. These are eliminated by integrating over the en- 
ergies of particles and holes generated in the Fermi seas 
during the process, taking into account thermal occupa- 
tions by / or 1 — / factors, where / is the Fermi function. 
As a result, one arrives at rates W which are labelled by 
molecular degrees of freedom only. (Of course, the rates 
still include a label for the participating junction(s).) 

To lowest order in Ht, one obtains sequential tunneling 
processes, which transfer one electron from a lead onto 
the molecule or vice versa. For example, the rate for 
a process changing the charge number of the molecule 
from n = to n = 1 by tunneling across junction a, and 
simultaneously changing the number of excited phonons 
from q to q' , is given by 

W°^,., = s{0, 1)^ \M,,,f + W - g]M- (7) 

Here, /a(e) = /(e — Ma) denotes the Fermi function for 
lead a. The symbol Mggi denotes the FC matrix ele- 
ments for a phonon transition q q', given by the over- 
lap of two harmonic oscillator wavefunctions (f>, spatially 
displaced by the distance Ax = \f2\losc- Here, i?osc = 
-\/ (a;2)q=o = {h/MujY/"^ defines the oscillator length. In- 
troducing q = min{gi, q^} and Q = maxjg'i, 52}, the ma- 
trix elements read 

/oo 
dx(t>l^{x + Ax)<t>q^{x) (8) 

= [sgn(g2 - q^T-"^ X^-'^e-^' ^LQ-''{X% 

where L^(a;) denotes the generalized Laguerre polyno- 
mial. 
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Our choice of jointly treating the spin-up and the spin- 
down states as the singly-occupied state n = 1 requires 
the inclusion of a spin factor s. The definitions s(0, 1) = 2 
and s(l,0) = 1 reflect the fact that rates for processes 
n = ^ 1 are twice as large as rates for n = 1 — + due 
to the spin-degeneracy of the state n = 1. The rate for 
the analogous process n = — > 1 can be obtained from 
Eq. (7) by swapping the spin factor and substituting fa 
by (1 - fa), i-e. 

W^^,,a = s{l,0)^ \M,,,f {1 - faied ~ W - q]hio)} . 

(9) 
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FIG. 3: Relevant contributions to the cotunneling rates in 
the U oo limit. For both charge states n = and 1 there 
are two incoherent contributions with either a spin-up or a 
spin-down electron being transferred across the molecule. 

Here, the spin factors are s(0, 0) = s(l, 1) = 2, which 
follows from an analysis of the relevant contributions as 
illustrated in Fig. 3. For cotunneling through the empty 
molecule, n = — > 0, there is an incoherent addition 
of the processes transferring a spin-up or a spin-down 
electron, respectively. For n = 1 ^ 1 the cotunneling 
transition may either leave the molecule's spin state in- 
variant or cause a spin flip, again resulting in two inco- 
herent processes. In both cases n = and 1, the rates 
of the two contributions are identical and hence can be 
absorbed into the spin factor. 

It is crucial to note that, in general, the cotunnel- 
ing rates cannot be directly evaluated from Eqs. (10) 
and (11). These expressions diverge due to second-order 



In next-to-leading order in Ht, cotunneling processes 
are generated. These transfer one electron from lead a 
to lead b, while the electronic occupation of the molecule 
changes only virtually in the intermediate state. Based 
on Eq. (6), the corresponding rates are obtained as 



(10) 



(11) 
I 

poles from the energy denominators. This divergence, 
related to the infinite lifetime of the virtual interme- 
diate state within a purely perturbative approach, has 
been pointed out before. ^^'^^ A regularization scheme has 
been developed in order to extract the correct cotunnel- 
ing rates, ^^'^"^ which we summarize here, deferring the 
details to Appendix A. 

The basis for the regularization of cotunneling rates is 
given by two observations, (i) Second-order perturbation 
theory clearly misses the fact that the intermediate state 
obtains a finite width ^ F due to the tunneling. This 
width should enter the energy denominator as an imagi- 
nary part, shifting the pole away from the real axis, (ii) 
For a specific cotunneling transition | i ) — > | / ) at finite 
temperature, the final state | / ) can also be reached from 
the initial state | i ) by two sequential processes. 

Specifically, the regularization proceeds as follows. 
First, a level width 7 ~ F is introduced in the energy 
denominators. With the poles shifted away from the real 
axis, the integral can indeed be carried out. The result- 
ing expression is then expanded in powers of the level 
width 7. The leading order term is proportional to I/7. 
When combining this with the prefactor FaFf, of the rates 
Eqs. (10) and (11), the overall order of this term is found 
to be identical to a sequential-tunneling term. Following 
point (ii), this term must be disregarded to avoid double- 
counting of sequential processes. The next-to-leading or- 
der in the 7-expansion is a 7*^ term. This term gives the 
regularized expression for the cotunneling rate. The full 
regularization procedure is carried out in detail in Ap- 
pendix A, and the resulting rates are given in Eqs. (Al) 
and (A3). 



E 



E 



Tf ed-e+{q' - q")hio 



faie) [l-/fc(e+[g-g']M] 



faie) [l~ftie+[q-q']hu;)] 



6 



It is interesting to note that these "cotunnehng rates" 
may become negative. This result can in fact bo ex- 
pected by considering the simple case of a single-level 
resonant tunneling model without charging energy. For 
this model, the exact result for the current in the T ^ 
limit is 



I = sgn(ey) 



2e TlTr 



E 



a arctan 



a\eV\-2ed 



(12) 

where we identify the index L (R) with a = +1 {a = —1). 
Expanding this in powers of F, we find that the second- 
order term corresponding to cotunneling is given by 



-^cot — 



Ir F l^^l 

h ^ "^(6^/2)2- 



(13) 



This expression indeed becomes negative, whenever the 
level is located inside the bias window, i.e. cotunneling 
reduces the sequential tunneling. Our regularization pro- 
cedure combined with the rate-equations approach ex- 
actly reproduces these expressions, as shown in detail in 
Appendix B. It is worth noting that, alternatively, the 
formalism developed by Konig, Schoeller and Schon may 
be used to circumvent the necessity of regularization, see 
e.g. Ref. 31. We find that additional corrections due to 
level shifts and broadening captured by their approach 
are irrelevant in the case F <C A;bT considered here. In 
this limit, both approaches are expected to give identical 
results. 

In the FC blockade regime, the dominant contribu- 
tions to the cotunnehng rates [Eqs. (10), (11), and (Al), 
(A3)] arise from highly excited virtual phonon states. Ac- 
cording to the Franck-Condon principle, the wavefunc- 
tion overlap reaches a maximum for transitions that oc- 
cur (approximately) vertically in the diagram Fig. 1(b). 
For transitions originating in the vibrational ground state 
g = 0, the "optimal" overlap is obtained for virtual states 
q" of the order of q" « A^. Deep inside the FC block- 
ade, the resulting poles e « ± X'^huj, see Eqs. (10), 
(11), are typically located well beyond the bias window. 
As a consequence of that, negative cotunneling rates are 
essentially irrelevant in the FC blockade regime. 

For an understanding of cotunneling rates it is im- 
portant to note that there is not a single virtual state 
q" resulting in maximal overlap, but rather an extended 



range Aq ~ A of excited states centered around q" 
see Fig. 1(b). By applying Sterling's formula to 



expression (8) 
obtain [Mqla^ 



for the FC 

|4 



matrix element Mn 



the 
we 



2j I « l/27rA2, so that the contribution 
from the intermediate state q" = [A^fiwJ is given by 

r^Ft, keT+eV 

A, and subsequent squaring 



the effective range Aq 
leads to an additional factor of A^, resulting in the ap- 
proximation 



TTh (ed + A2fiw)2 



(14) 



valid for < eV <Si }?hw. A more detailed derivation of 
this expression is given in Appendix C. By contrast to 
the sequential-tunneling rates, Eq. (14) establishes that 
cotunneling rates are not exponentially but rather alge- 
braically suppressed. In particular, for = 0, we have 



B. Current and shot noise 

As shown, e.g., in Refs. 10 and 32, the steady-state 
current and the current noise can be computed via the 
rate-equations formalism. In the first step, the steady- 
state occupation probabilities are obtained from the 
rate equations 



0= — ^ 
dt 



q'q 



^ q " qq' 



(15) 



The last term is added on a phenomenological basis, 
and describes relaxation of the vibrations towards the 
equilibrium distribution P^'i = e''^'^/''^'^ {1 - e-'^^/'^f^) 
on a time scale r. In the following, we will mainly in- 
vestigate the two limiting cases of equilibrated phonons 
(r ^ 0, distribution fixed to equilibrium), and unequi- 
librated phonons (r ^ oo, nonequilibrium distribution 
entirely determined through tunneling dynamics). 

It is convenient to reexpress Eq. (15) in terms of a 
matrix equation of the form = WP, where W is a co- 
efficient matrix containing all rates, and the vector P 
consists of all steady-state probabilities. [See Appendix 
D for further details.] For a systematic expansion in or- 
ders of the tinincling, we write W = W^^' + W^^) + . . . as 

well as P = P(") + P(^) H . Substituting this into the 

rate equation, we obtain up to second order the equations 



W(i)p(0) = 0, 
wWp(i) =-W(2)p(0). 



(16) 
(17) 



The normalization condition for P leads to trP(°) = 
^ . P^°^ = 1 and tr P(i) = 0. Here and in the follow- 
ing, the "trace" of a vector always denotes the sum of 
its components. Similarly, we expand the steady-state 
current 

/ = /(!) +7(2) + ... (18) 

= tr W«P(0) + tr [W(^)p(i) + Wf P(0)] + • • • 

Here, W/ denotes the reduced coefficient matrix contain- 
ing only current-carrying processes, see Appendix D. The 
lowest-order current /*^^^ is the sequential-tunneling cur- 
rent. The second-order contribution I^"^^ contains the co- 
tunneling current and corrections to the sequential cur- 
rent due to changes of the probability distribution by 
cotunneling. 
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A method for evaluating the current shot noise 

/oo 
dTe--[(7(T)7(0))-(7)2] (19) 
-co 

within the rate-equations formaUsm has been developed 
by Korotkov in Ref. 32. This method can be general- 
ized to include higher-order processes beyond sequential 
tunneling, as we show in detail in Appendix D. It is 
an important question whether the inclusion of negative 
"cotunneling rates" into this noise formalism is justified. 
In Appendix B, we explicitly verify for the resonant- 
tunneling model that this is indeed the case. The exact 
solution for the resonant-tunneling model allows for an 
evaluation of the zero-frequency noise via the relation 



S=^ Jdsi^^TiEJUE) [l-fa{E)] 



(20) 



+ T{E) [1-T{E)] [ME) -ME)] 



see Refs. 18,33. We find that the leading and next-to- 
leading order expansion in r/fc^T of the resulting noise 

is exactly reproduced by the generalized Korotkov for- 
malism including negative cotunneling corrections. 



V. CURRENT- VOLTAGE CHARACTERISTICS 

A. Linear Conductance 

For discrete electronic levels weakly coupled to two 
electrodes via tunneling junctions, the linear conduc- 
tance is known to develop peaks whenever an electronic 
level is aligned with the Fermi energies of the leads. 
These peaks are well described by sequential-tunneling 
processes within the rate-equations formalism. Away 
from the peaks, the sequential-tunneling conductance 
is exponentially suppressed in the parameter ed/ksT, 
and cotunneling contributions become important. In 
the following, we show that strong electron-phonon cou- 
pling destroys this gate-controlled crossover between se- 
quential tunneling and cotunneling at low temperatures 
ksT -C hu, rendering cotunneling dominant over wide 
parameter ranges. 



1. Sequential-tunneling conductance 

We start by considering the lowest-order contributions 
to the conductance, in the limit of low temperatures F <C 
fcsT ^ huj. In this case, only the vibrational ground 
state is occupied, and any real phonon excitations are 
negligible. This only leaves the rates Wqq and Wqq for 
consideration. Employing the rate-equations approach to 
this simple situation, the linear conductance for infinite 
U is obtained as 



^seq — 



2e2 FrF 



Li 77, 



fi^d) 



h TL+TRl + }{ed) 



(21) 



In the absence of vibrations, this result can be found, 
e.g., in Ref. 34. The crucial point here is the fact that the 
presence of vibrations leads to an exponential suppression 
of the linear conductance in A^, even at the peak value. 
As discussed in Ref. 34, it is interesting to note that 
in distinction to the U = Q case, the conductance peak 
position is not at = but varies with temperature as 
= A;BTln2/2. This shift originates from the fact that 
the sequential-tunneling rates n ^ 1 arc twice as 
large as the rates n = 1 ^ due to the spin degeneracy. 



2. Cotunneling conductance 

We now consider the next-to-leading order contribu- 
tions to the conductance. Remaining in the limit of low 
temperatures T <^ ksT <^ huj, the only two additional 
processes are elastic cotimneling through the empty and 
the singly-occupied molecule with corresponding rates 
^oo-aa' ^^'^ ^oo-aa'^ respectively. The resulting conduc- 
tance due to elastic cotunneling is given by 



d 



dV 

2e^ FrT 



1=0,1 



nn 

00;LR 



117"" 



(22) 



v=o 



Li- R 



h l + f{ed) 



1 - fi^d) 

{ea + A2/ia;)2 ^ (e^ - X'^hu;)^ 



2/(ed) 



Here, the apparent poles at Cd = ±X^hco are exponen- 
tially suppressed due to the Fermi function factors. In 
fact, they are an irrelevant artefact of our simple ap- 
proximation in Eq. (14). The correct treatment via the 
regularization procedure eliminates this flaw. 

By contrast to the ^ e~^^ suppression of the 
sequential-tunneling conductance, the cotunneling con- 
ductance is only algebraically suppressed as ~ A~^. 
This is a strong indication that, in the regime of strong 
electron-phonon coupling, cotunneling will dominate the 
linear conductance, even at the conductance peak. Rep- 
resentative numerical results confirming this expectation 
are shown in Fig. 4(a). For ksT <C fkv, the cotunnel- 
ing conductance Gcot is dominant. In this interesting 
case, the conductance develops a kink feature at = 
broadened over a gate- voltage range given by ksT, which 
can be understood from Eq. (22). Away from the peak 
maximum, Gcot roughly decays as ~ (e^ -|- const. )~^, see 
Fig. 4(b). This is consistent with Eq. (22). For larger 
temperatures ksT > 0.2hoj, real excitation of phonons 
becomes important. Such excitations, which we have ne- 
glected in our approximations (21) and (22), open up 
additional channels for sequential tunneling. The in- 
crease of sequential contributions is evident from the 
high- temperature curve in Fig. 4 (in red/dark grey color), 
which shows the recurrence of a sequential-tunneling 
peak on top of the cotunneling background. 
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w/cotunneling, ksT = O.ihij 

w/ cotunneling, fesT = O.SRtj 

w/cotunneling, ksT=0.5F^uj 

w/o cotunneling, ksT = O.LHuj 

w/o cotunneling, ksT = O.Shij 
w/o cotunneling, kBT = 0.5tiLLi 



FIG. 4: (Color online) Linear conductance as a function of 
gate voltage in the FC blockade regime (A = 4, Fa = 0.02/icj). 
(a) A comparison of the sequential-tunneling conductance 
(dashed lines) and the conductance including cotunneling 
corrections (solid lines) clearly shows the dominance of co- 
tunneling contributions at low temperatures. At e^j « 0, 
cotunneling results in a conductance kink for fesT <C huj. 
Sequential contributions become relevant for temperatures 
fcsT > 0.2hui that allow thermal-induced occupation of ex- 
cited phonon states, (b) The plot of G~^^^ reveals the in- 
verse quadratic decrease of the cotunneling conductance as a 
function of gate voltage. 



B. Nonlinear current- voltage characteristics 

1. Sequential current 

For strong electron-phonon coupling, we have estab- 
lished in Ref. 11 that the nonlinear IV develops a low- 
bias gap, essentially eliminating the conventional crossing 
at the degeneracy point, cf. Fig. 7(b). For equilibrated 
phonons, it is not difficult to estimate the sequential- 
tunneling current in the FC blockade regime. The main 
input is that (i) for low temperatures fc^T <C huj and 
strong relaxation, the initial vibrational state is the 
phonon ground state q — for each tunneling pro- 
cesses, and (ii) the squared Franck-Condon matrix ele- 
ments I Moq I strongly increase with q up to values g ~ . 

Due to (ii), the current is dominated by tunneling 
events with maximal change in the phonon number, 
Aq = [\eV\ /2hLu\, as long as \eV\ /2hLu < A^. Hence, 
(i) allows us to approximate the current via the rates 



/.eq«Sgn(l/)-- 



6 2F(jr(j' 



Aq\ 



(23) 



with (a, a') = (i, R) or (i?, L) depending on the bias sign. 
Here, we have assumed = for simplicity, which can in 
principle always be arranged by adjusting the gate volt- 
age. After averaging over steps, one obtains for voltages 
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eV [hLo) 
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FIG. 5: (Color online) Comparison between the sequential- 
tunneling current [black curve] and the asymptotic formula 
(24) [red/grey curve] for equilibrated phonons in the Franck- 
Condon blockade regime, (e^ = 0, ksT = 0.05fta;, Fl = F^) 



larger than temperature 



-X^ •^\eV\/hu>- 



h 2Ta + Ta' T{\eV\ /ihuj + 1/2) 



sgn(T/), (24) 



where the denominator involves the Gamma function 
F(x). 

The comparison between the asymptotic form Eq. (24) 
and the actual rate-equations solution is depicted in 
Fig. 5. Deviations of the asymptotic form from the actual 
solution are observed when the inequality 2\eV\ <^ X^hjj 
gradually breaks down. Then, the true current assumes 
higher values as compared to Eq. (24) , since assumption 
(ii) is no longer valid, i.e. not only one but several phonon 
channels yield significant contributions to the current. 

For unequilibrated phonons, where FC blockade leads 
to self-similar avalanches of electrons, ^^'^^ an estimate 
for the sequential current is less trivial. The avalanche 
dynamics makes it necessary to include several phonon 
channels. In particular, following the arguments of 
Ref. 17, the mean electron-number per avalanche as 
well as the mean number of sub-avalanches within each 
avalanche will depend on bias and gate voltage, as well 
as on the coupling strength A. 



2. Nonlinear IV including cotunneling corrections 

We now investigate the cotunneling corrections to 
the current-voltage characteristics in the FC blockade 
regime. Representative results for both equilibrated and 
unequilibrated vibrations are depicted in Figs. 6 and 7. 

We first focus on the case of equilibrated phonons. The 
IV plots for ed — 0, Fig. 6(a) and (c), illustrate the 
importance of cotunneling contributions well beyond the 
linear response regime. Cotunneling corrections signifi- 
cantly increase the current for bias voltages \eV\ < X^hui, 
the reason being again the increased overlap for oscillator 
states in cotunneling processes with highly excited inter- 
mediate states. It is a crucial point, however, that cotun- 
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eV [fk^) eV (fiw) 

FIG. 6: (Color online) Nonlinear current-voltage characteris- 
tics in the FC blockade regime. Each graph shows the steady- 
state current j'^-' + Z*-^' (solid curves), and the sequential 
current J*-^^ (dashed curves) for comparison. The electron- 
phonon coupling and vibrational relaxation strength are cho- 
sen as (a) A = 4, equilibrated, (b) A = 4, unequilibrated, (c) 
A = 5, equilibrated, (d) A = 5, unequilibrated. The results 
unambiguously indicate the dominance of cotunneling contri- 
butions in the low-bias FC suppressed regime. (Parameters: 
€d = 0, ksT = Omhcj, Fa = 0.02fi,cj) 



neling does not lead to a breakdown of the FC blockade. 
The blockade remains intact, since cotunneling rates are 
still algebraically suppressed at strong electron-phonon 
coupling. This suppression remains strong due to the 
rather large exponents in the scaling with A, see Eq. (14). 
An additional change caused by the dominance of cotun- 
neling rates is a drastic reduction of step features in the 
current, see especially Fig. 6(c). When plotting the dif- 
ferential conductance dl/dV on a logarithmic scale, see 
Fig. 7, one observes that small steps persist on top of 
the (much larger) cotunneling background and form the 
usual system of vibrational sidebands. 

An important low-bias feature beyond the sequential- 
tunneling approximation consists of two steps in dl/dV 
for eV = zthu!. These steps mark the onset of inelas- 
tic cotunneling. The emergence of such steps is not re- 
stricted to the FC blockade, and has been discussed, e.g., 
in Ref. 8. Additional steps at bias voltages equal to larger 
integer multiples of hio can be observed when going to 
even lower temperatures than chosen in Figs. 6 and 7. 

For unequilibrated phonons, we find that cotunneling 
corrections again dominate the FC blockade region of the 
current-voltage characteristics, see Fig. 6(b), (d). How- 
ever, relative to the equilibrated case, the corrections are 
less drastic. In particular, vibrational steps in IV remain 
well visible, and the bias range where cotunneling domi- 
nates is narrower than in the equilibrated case. This can 
be explained by the persistence of avalanche dynamics 
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FIG. 7: (Color online) (a) Differential conductance dl/dV as 
a function of gate and bias voltage for A = 4 and Ta = 0.02/iaj, 
plotted with a logarithmic color scale. The low temperature 
plots (left-hand side) show a step in dl /dV due to the onset 
of cotunneling (horizontal feature at eV — ±huj). For unequi- 
librated phonons, inelastic cotunneling gives rise to phonon 
absorption-induced vibrational sidebands inside the Coulomb- 
blockade region (whose borders are marked by dashed lines). 
For comparison, the plots in (b) show the corresponding re- 
sults within the sequential-tunneling approximation which 
misses the absorption-induced sidebands. (G' — xF^^n^) 



generated by sequential tunneling. Within the FC block- 
ade, important contributions to the current stem from 
rather featureless cotunneling events. However, phases 
of pure cotunneling are still interrupted by avalanches. 
During those, the system attains phonon excitations high 
enough to make sequential tunneling dominate over co- 
tunneling. Hence, sequential tunneling retains a signifi- 
cant role in the unequilibrated case. 

Fig. 7 demonstrates an additional striking consequence 
of cotunneling: For unequilibrated phonons, additional 
vibrational sidebands appear in the Coulomb-blockade 
region. The mechanism which generates these addi- 
tional sidebands is the following. The Coulomb-blockade 
regime, where transport is conventionally dominated by 
cotunneling, allows for excitation of vibrations through 
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inelastic cotunneling as soon as the bias voltage ex- 
ceeds the phonon energy. For weak vibrational relax- 
ation, such excitations decay slowly and therefore re- 
main relevant for subsequent tunneling events. The cru- 
cial point is that the usual exponential suppression for 
sequential-tunneling in the Coulomb blockade applies to 
processes which leave the vibrational state unchanged,^^ 
but not to processes gaining energy through phonon de- 
excitation. Thus, phonon excitations (induced by in- 
elastic cotunneling) facilitate the recurrence of sequen- 
tial tunneling accompanied by phonon deexcitation in 
the Coulomb-blockaded region. Furthermore, the onset 
of inelastic cotunneling at \eV\ = hw and the resulting 
absence of vibrational excitations for lower bias voltages 
explain why absorption-induced sidebands must termi- 
nate at \eV\ = hw and cannot enter the low-bias region. 

This interplay between inelastic cotunneling and se- 
quential tunneling results in the observed absorption- 
induced vibrational sidebands,^" resembling experimen- 
tal findings in transport through suspended carbon- 
nanotubes.""^^ We emphasize that the absorption- induced 
vibrational sidebands are not specific to the case of strong 
electron-phonon coupling, and point out that their pres- 
ence may be exploited to estimate the vibrational relax- 
ation rate in molecular junctions. 

Interestingly, phonon excitation due to inelastic co- 
tunneling is the only viable mechanism for absorption- 
induced vibrational sidebands up to second order in the 
tunneling. It is crucial to note that thermally induced 
phonon excitations will in general not lead to visible ad- 
ditional sidebands, but rather to a smooth smearing of 
the Coulomb diamond, see Fig. 7(b). This is explained 
by the fact that vibrational sidebands are clearly visible 
only for /c^T <C huj, and in this case thermally activated 
phonon excitations are negligible. 



VI. CURRENT SHOT NOISE 

The current shot noise serves as an important tool for 
acquiring information about the transport dynamics of 
the system beyond the steady-state current. In partic- 
ular, we are interested in the Fano factor F and noise 
spectrum S{u!) for unequilibrated vibrations to determine 
the fate of the self-similar avalanches in the presence of 
cotunneling corrections. We first turn to the regime of 
strong vibrational relaxation and show that even in this 
case, cotunneling leads to interesting new effects. 



A. Zero-frequency noise for equilibrated 
vibrations: Telegraph noise 

As shown in Ref. 11, sequential tunneling results in 
conventional sub-Poissonian Fano factors for equilibrated 
phonons, see the dashed curves in Fig. 8. Electrons are 
transferred one by one, interacting with each other only 
via the charging energy J7, which leads to the suppres- 
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FIG. 8: (Color online) Excess noise Fano factor for equili- 
brated phonons with A = 4 (top row) and A = 5 (bottom row). 
For nonzero ed we observe large super-Poissonian Fano factors 
due to current telegraph noise. {ksT — 0.05huj, Ta ~ 0.02huj) 



sion of F below 1. As depicted in Fig. 8 by the black 
curves, the inclusion of cotunneling processes radically 
changes this picture. The overall Fano factor remains 
sub-Poissonian for = 0. Remarkably, it switches to 
super-Poissonian behavior for nonzero in a certain bias 
range. The magnitude of F strongly depends on the cou- 
pling strength A, as well as gate and bias voltage, and 
reaches values as large as 100 for A = 5. 

More insight into this surprising enhancement of noise, 
which even occurs in the low-bias regime where no 
phonons can be excited, is obtained by Monte-Carlo sim- 
ulations. For the parameters of the large maximum of 
the Fano factor in Fig. 8 (A = 5, = O.bhuj), Fig. 2 
shows the coarse-grained current as a function of time. 
These data clearly show the slow switching between a low 
and a high current state, on top of much more rapid dy- 
namics within each current state. The underlying reason 
for this telegraph noise is the role inversion of sequential- 
tunneling and cotunneling rates in the FC blockade, i.e. 
rates for cotunneling are larger than those for sequential 
tunneling. 

In the following, we give a detailed explanation of 
the telegraph behavior. Denoting the molecular state 
by I n, g ) , a straightforward inspection shows that for 
ed = O.bhuj and eV = O.Shuj, the only relevant transitions 
are (1) I 0, ) ^ 1 1, ) and I 1, ) ^ I 0, ) due to sequen- 
tial tunnehng, and (2) | 0, ) ^ | 0, ) and 1 1, 1 ) ^ 1 1, 1 ) 
due to cotunneling. As explained in Section IV, the co- 
tunneling transitions are much faster than the sequen- 
tial transitions in the FC blockade regime, such that 
Woo, ^00 > ^00 > Woo- The key point in the emer- 
gence of the telegraph behavior is the fact that the rates 
for cotunneling through the empty and through the singly 
occupied dot are generally different, Wqq ^ IFqo^- This 
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can be understood from Eqs. (10) and (11), noting that 
Cd enters the energy denominator with a positive sign for 
n = 1, but with a negative sign for n = 0. For = 0, 
the two cotunnehng rates are identical, and indeed no 
supcr-Poissonian noise is detected in this case. However, 
for nonzero ea, the two rates become different. As a con- 
sequence, the rapid cotunnehng dynamics in the charge 
states n = 0, 1 leads to different mean currents eWgo and 
eW^Q , respectively. The rare sequential transitions effect 
the slow switching between these two current states. 

The Fano factor resulting from the current telegraph 
noise is directly related to the excess charge transmitted 
in the high-current state during the typical residence time 
in the corresponding charge state, i.e. F ~ AWcot/Wseq. 
Here, AWcot = |l^oo — W^oo^l denotes the difference of 
the relevant cotunnehng rates. Due to the role inver- 
sion of cotunnehng and sequential tunneling, we have 
Wcot ^ W^seq SO that the resulting Fano factor may be- 
come large. A more quantitative calculation of the Fano 
factor can be achieved by employing the generalization 
of Korotkov's method'^^ described in Appendix D. Along 
these lines, we can derive approximate expressions for the 
zero-frequency shot noise and its Fano factor, 



S{f = 0) « 4e^ 



F « 6 



(25) 
(26) 



see Appendix E for details. This approximation is valid 
whenever Wo"o°, W^^ > W^^, W^^ and W^^ ^ W^^. 
Counting the orders of cotunnehng versus sequential- 
tunneling rates in the numerator and denominator of the 
Fano factor, we find that the Fano factor is indeed of the 
order of AVFcot/Wscq- 



B. Zero- frequency noise for unequilibrated 
vibrations: Persistence of avalanches 

For unequilibrated vibrations, the sequential-tunneling 
dynamics in the FC blockade regime is described by 
self-similar avalanches of electrons. ^^'^^ The large Fano 
factors due to sequential tunneling shown in Fig. 9 
(dashed curves) reflect the large number of electrons per 
avalanche. In the = case, the Fano factor is ex- 
actly given by = 2{Ni), i.e. twice the mean number 
of electrons per level-0 avalanche. More generally, the 
Fano factor can be obtained by evaluating the fluctua- 
tions of the waiting times ti between level-0 avalanches, 
and the fluctuations of the electron number Ni in each 
avalanche, 



F = (N.,) 



{tl)-{U? , {Nf)-{N,)' 



(27) 



Before turning to the cotunneling-induced corrections, 
we study the effect of nonzero e^, which was not discussed 
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FIG. 9: (Color online) Excess noise Fano factor for unequi- 
librated phonons with A = 4 (top row) and A = 5 (bottom 
row) with ksT = 0.05/ia; and Ta ~ 0.02huj. Due to the co- 
existence of cotunnehng and avalanches, the Fano factor is 
reduced as compared to the sequential-tunneling approxima- 
tion. The super-Poissonian peak observed at low bias for 
A = 5 and ed 7^ signals the occurrence of current telegraph 
noise. 



in Refs. 11 and 17 in detail. Interestingly, a nonzero ed 
can lead to even larger Fano factors as compared to the 
Cd = case, see Fig. 9. Especially large Fano factors 
are observed for low biases above the threshold eV > hui 
whenever the electronic level is roughly aligned with 
the Fermi energy of the right or left lead. 

At first view, the resulting Fano factor enhancement 
is surprising since the mean number of electrons per 
avalanche is decreased for such a level configuration. For 
an understanding of this decrease, consider the case of 
level alignment with the left Fermi energy. Then, phonon 
excitations can only occur through tunneling in the right 
junction. Events in the junction L can only decrease 
the phonon state (or leave it unchanged). Thus, every 
second tunneling event will inhibit phonon excitations. 
With avalanches benefitting from the possibility of reach- 
ing highly excited phonon states, the mean number of 
electrons per avalanche is reduced as compared to the 
case where the level is centered between the two Fermi 
levels. This decrease in (Ni) is confirmed by Monte-Carlo 
simulations. 

The somewhat counterintuitive result that the Fano 
factor nevertheless becomes larger can be traced back to 
the significantly enhanced fluctuations in waiting times 
ti, which directly contribute to the Fano factor, see 
Eq. (27). In the following, we explain this enhancement 
of waiting-time fluctuations. Any generation-0 avalanche 
may terminate either in the charge state n = or n = 1, 
and the corresponding probabilities are denoted po and 
Pi. After the avalanche termination, the system un- 
dergoes a waiting period ti. The crucial point is that 
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for level alignment with one Fermi energy, this waiting 
time strongly depends on the charge state n. Starting 
in the state n = 0, transitions in the left junction do 
not allow phonon excitations. By contrast, for the ini- 
tial state n = 1 and transitions across the right junc- 
tion, the excitation of phonons is favored (due to the 
increase of FC matrix elements towards highly excited 
vibrational states). Accordingly, the relevant two rates 
after avalanche termination are dramatically different, 
Wq^ <^ Wqi, and the waiting times in the two cases 



will differ accordingly, 



(28) 



Consequently, the probability distribution of waiting 
times is not given by a single exponential anymore, but 
is well approximated by 

piU) ^ PoW°^ eM-KoU] +pi<i°exp[-<°t,]. (29) 

While the ed = case (with a simple exponential distri- 
bution) results in — (ij)^) / (ti)^ — 1, we now find 



Pq 



Po 



(30) 



vahd for po '3> Wqo /^oi*- ^y Monte-Carlo simulations 
we confirm that po satisfies 1 > po ^ ^^00/^^0° for re- 
alistic parameters. Therefore, Eq. (30) explains the en- 
hancement of Fano factor due to the increased waiting- 
time fluctuations. 

The inclusion of cotunncling contributions to the cur- 
rent shot noise results in the black curves depicted in 
Fig. 9. Again, at low biases which do not allow real 
phonon excitations, we observe traces of telegraph noise. 
[See the A = 5 plots for nonzero in Fig. 9.] For higher 
biases which facilitate avalanche dynamics, we find that 
cotunneling generally reduces the Fano factor as com- 
pared to pure sequential tunneling. The reason for this 
reduction is the coexistence of electron avalanches (lead- 
ing to large Fano factors) with cotunneling, which re- 
sults in Fano factors of the order of unity away from 
the telegraph noise regime. We emphasize that despite 
the cotunneling-induced noise reduction, the overall Fano 
factor remains large compared to unity, signalling the 
persistence of electron bunching in avalanches. 



C. Noise spectrum for unequilibrated vibrations: 
Self-similarity of avalanches 

Probing the self-similarity of avalanches requires one 
to go beyond the zero-frequency noise. In Refs. 11,17 
two of us have shown that the characteristic signature 
of the self-similarity within the sequential-tunneling ap- 
proximation is an approximate power-law behavior of the 
noise as a function of frequency f,S^ f^^^'^- In Figure 
10 we present the results for the noise power spectrum 
including cotunneling corrections. 
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FIG. 10: (Color online) Noise power spectrum for unequili- 
brated phonons with A = 4 (left) and A = 5 (right) with 
ksT = Omhuj, Fa = 0.02huj, €d = 0, and eV = 3hio. The ap- 
proximate power-law behavior reflecting the self-similarity of 
avalanches persists in the presence of cotunneling. The domi- 
nance of cotunneling during sequential waiting times provides 
a long-time (low-frequency) cutoff for the power law. 



The central result is that the power-law scaling remains 
valid. Only the frequency range over which it applies 
is reduced by cotunneling. In particular, we observe a 
larger low-frequency cutoff when taking cotunneling into 
account. This result is consistent with the picture of co- 
existence of cotunneling and sequential avalanches: In 
the sequential-tunneling picture, the onset of the power 
law at low frequencies is essentially given by the strongly 
suppressed rate for escaping from the vibrational ground 
state, i.e. / WqI^^^. Going beyond sequential tunnel- 
ing, we have shown that transitions close to the vibra- 
tional ground state are dominated by cotunneling which 
does not contribute to the avalanche dynamics. Accord- 
ingly, the onset of self-similarity and hence the power-law 
behavior is shifted towards higher phonon levels with cor- 
responding higher rates. This is reflected by a larger low- 
frequency cutoff for the power law in the noise spectrum. 



VII. CONCLUSIONS 

In conclusion, we have developed a complete theory of 
the Franck-Condon blockade including cotunneling cor- 
rections. This regime, caused by a strong coupling to the 
bosonic degrees of freedom of molecular vibrations, em- 
phasizes the stark contrast between transport through 
molecular junctions and more conventional nanostruc- 
tures such as quantum dots. In particular, the intriguing 
finding that electron-phonon coupling may result in a 
significant bunching of electrons into a hierarchy of self- 
similar avalanches dominating the transport, constitutes 
a novel collective effect unknown from electron transfer 
through quantum dots. 

In distinction to the Coulomb blockade scenario in 
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quantum dots, where the crossover between sequential 
tunnehng at the conductance peaks and cotunnehng 
away from the peaks may be probed by tuning the gate 
voltage, we find that the FC blockade regime features 
the dominance of cotunneling contributions in the entire 
low-bias region. The underlying reason for this result is 
the enhanced overlap between vibrational state, which 
is exploited by highly excited vibrational states in co- 
tunneling processes. We have shown that the FC block- 
ade and avalanche-type transport persist despite the co- 
tunneling corrections. Furthermore, we have established 
that cotunneling leads to a number of remarkable effects 
missed by the sequential-tunneling approximation. In 
summary, our central results are: (i) Apart from horizon- 
tal dl/dV steps due to inelastic cotunneling, the most im- 
portant consequence of cotunneling is the appearance of 
additional absorption-induced vibrational sidebands in- 
side the Coulomb-blockade regime, relevant for unequili- 
brated vibrations. This may also serve as an experimen- 
tal measure of the vibrational relaxation rate, (ii) The 
role inversion between cotunneling and sequential tun- 
neling may cause current telegraph noise whenever the 
electronic level is not centered between the left and right 
Fermi energy. This is reflected in strongly enhanced Fano 
factors, independet of the vibrational relaxation strength. 

For a more detailed summary of results we refer the 
reader to the summary provided in Section I. In clos- 
ing, we emphasize the appeal of realizing shot noise 



measurements in experiments with single-molecule de- 
vices. These could provide essential information about 
the transport mode, which cannot be accessed by current- 
voltage characteristics. Experimental efforts in this di- 
rection are under way in several research groups. 
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APPENDIX A: REGULARIZED COTUNNELING 
RATES 

The regularized cotunneling rates can be expressed in 
terms of digamma and trigamma functions, 



0^0 _ ^-^a-^b 

27r?i 



q^q ,ao 



^\Mq,rMqrf J {^a, Hb - [q- q']hv,ed- [q-r]hu)) 

r 

+ Mq,rM*^M*,^MqJ {na, lib-[q- q']t^, id - [q - r]fiw, ed-[q- s]t>io) 



(Al) 



r^s 



q^q ,0,0 



2Trh 



J2 \Mq,rMqr\^ J{iia,Hb - [q - q']Jv^,ed + [q' - r]huj) (A2) 

r 

+ ^ Mg,rM*^M*^MgsI {l^a, l^b-[q- q']fiu), ea + [q' - r]fuv, €4 + [q' - s]nco) 

Here, the regularized expressions / and J can formally be obtained by evaluating the 7^0 limit after subtracting 
1/7 terms. 



I{EuE2, ei, €2) = lim Re / dE f{E - Ei) [1 - f{E - E^)] 

7^0 J i/ — ei 

nB{E2 - El) 



1 



ei - £2 



?>/ _E - £2 + n 

Re <! V(l/2 + i(i[E2 - ei]/27r) - ^(1/2 - z/3[^2 - e2]/27r) 



(A3) 



- V(l/2 + il3[Ei - ei]/27r) + V(l/2 - i/3[£i - e2]/27r)| 



J{Ei,E2,e) = lim 
7— >o 



J dEf{E-Ei)[l-f{E-E2 



1 



{E - e)2 + 72 



- 0(1/7) 



(A4) 



^nB{E2 - El) lm{i,' (1/2 + i(3[E2 - e]/27r) - V' (1/2 + if3[Ei - e]/27r)} . 

ZTT 
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Our numerical implementation of the polygamma func- 
tion ip has been based on a routine written by Kolbig, 
available in the CERN Program Library. ^^'^^ 



APPENDIX B: EXACTLY SOLVABLE MODEL 

We consider the simple case of a single spin-degenerate, 
noninteracting level coupled to two leads. Due to the 
absence of interaction in this model, spin only results in 
trivial factors of 2 and may be neglected. We have picked 
this trivial model, since it is exactly solvable, e.g. using 
a one-particle scattering approach. 



Current 



The current is given by 

TLiE)TR{E) 



dE: 



A4E)[fL{E) - UE)], 
(Bl) 



TLiE)+TR{E) 
wfere denotes the spectral function of the dot, 

T{E) 



Ad{E) = -2lmG^{E) = 



[E-ea]^ + [T{E)/2]^ 



, (B2) 



and the energy broadening of the level is T{E) = rL{E) + 

TiiiE) with Ta{E) = 2TTpa{E) \ta\^. (Note that any pos- 
sible energy shift due to the tunnel coupling has been 
absorbed into e^.) In the wide-band limit and for sym- 
metric voltage splitting, the solution reads 



dE- 



HTl+TrJ ^"[S-e<i]2 + [r/2]2 
e TlTr 1 



ifUE) - fn{E)] 



(B3) 



Here, V denotes the digamma function. The exact result 
can be expanded in orders of F, which gives 



rti L + t R 



(B4) 



a ^ ^ 

For the weak-coupling case we can alternatively apply 
the rate-equations approach. In this case, the stationary 
current is calculated via I = II = —Ir with 



T _ 7-(seq) , 7-(cot) 



Here, (W^ ) gives the rate for a sequential-tunneling 
process in junction a where one electron enters (leaves) 



the dot. Similarly, W^r {Wp^"^ 



RLj 



denotes the rate for a 
cotunneling process with initial and final charge state n 
where one electron is coherently transferred from the left 
to the right lead (or vice versa). The stationary prob- 
abilities for the empty and occupied dot Pq and Pi are 
obtained from the rate equations 



: 



PiW- 



PoW^ 



Po + Pi^ 1, (B6) 



+ Wp- It is important to note that in 



dPo 
dt 

where 

this case, there are no cotunneling contributions to the 
rate equations: Elastic cotunneling does not change the 
charge state of the dot, and hence does not affect the 
probability distribution. 

Sequential-tunneling contributions. In the first step, 
we calculate the stationary probabilities and the 
sequential- tunneling current. The transition rates ob- 
tained from Fermi's golden rule are 



1 



From substituting this into Eq. (B6) one obtains 

w- _ E„ r„[i-/„(e,)] 



(B7) 
(B8) 

t L + t K 

and upon substituting into Eq. (B5), the sequential- 
tunneling current is found to be 

TlTr 



Pn = 



Pi = 



w- + W+ 
w+ 



Tl + Tr 



w- 



w+ 



/(seq) _ ^_ 



[/L(ed) - /fl(e<i)] 



(BIO) 



Indeed, this expression is indentical with the leading- 
order current contribution from the expansion of the ex- 
act solution, see Eq. (B4). 

Cotunneling contributions. — In the second step, we cal- 
culate the cotunneling current by evaluating the cotun- 
neling rates 



ah 





1 


1 de 





2TTh 

After regularization, we find 
13 TaYb 



/4e)[l-/fe(e)]. (Bll) 



WT 



47r2 h 
xIm{V''(i+i|;[/i6 



"s(M6 - Mo) 



4,' 



^ + i-^\u 



(B12) 

)}■ 



Substituting this together with Eqs. (B8), (B9) into the 
expression for the cotunneling current from Eq. (B5), and 
using nsix) + nB{—x) = — 1 leads to 



(B5) /(cot) =|r^r«Aim^a^' fl+^£.[e, + aeV/2]\ 

n. ^ ' 



(B13) 

By comparison with Eq. (B4) one verifies that this is 

identical with the next-to-leading order term of the ex- 
pansion of the exact solution. 
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2. Zero-frequency noise 

Our starting point for the exact noise calculation is the 
two-terminal formula for the zero-frequency noise, 

S=^J dE^^^TiE)UE)[^-ME)] (B14) 
+ T{E) [1-T{E)] [U{E)-fn{E)fy 

see Refs. 18,33. The transmission coefficient T{E) is 
known to be proportional to the spectral function Aci{E). 
The exact relation between them can be extracted easily 



by comparing Eq. (Bl) to the Landauer formula 

{I) = lJ dET{E) [h{E) - fnm . (B15) 

This leads to T{E) = Ad{E)TLTR/T. 

The integrations can be carried out by invoking the 
Fourier transform for the Lorentzian. The energy- 
integration then turns into a Fourier transform of the 
Fermi factors, which may be carried out by contour 
integration. The remaining time-integration leads to 
polygamma functions. The exact result for the zero- 
frequency noise is given by 



sJ-^[ ^^^^^g + coth(;jeF/2) Im (1 + z A [e, + aeV/2] ' 



An 



+ ^fff^^ ^ cotHPeV/2) Im 5^ a^' Q + z^ [e, + aeV/2] + ^) 



7rr3 



(B16) 

Expanding this expression in F, we find in leading and next-to- leading order 

^(,) = 2e^| 2^r^rH(ri + r|) ^^^^^^^ ^ ^^^^^^ _ ^^^^^^^^^^^^y^ (bi^^ 

^(2) ^ A coth(/3eF/2) Im^a V' + i^[e, + (B18) 



Again, wc can compare this expansion with the results 
from the rate-equations approach. Employing the for- 
malism described in Appendix D, we obtain the occupa- 
tion probabilities 



and the vectors 



Eara[l-/a] 



the steady-state current 



FrF 



J- L + J- _R 



R 



Tl+Tr 



Ul - Ir) + WlR - Wrl 



(B20) 



urr = 



URL 



YL = 



Yr 



Wl = 



Pn{^L.fL/f^ + WLR + Wrl) 
Pi{Tl[1 - .fL]/h + WLR + wrl) 

PoO^RfR/fl + WLR + Wrl) 
Pl(rfl[l - fR\/h + WLR + Wrl) 

_ ( Pq{wlr + Wrl) 
\Pi{wlr + wrl) 

-PiTl[1 - hyh + PoiWLR - WRL 
Po^LfL/h + PiiWLR-WRL) 



PiTrII - fR]/h + Pa{wLR - wrl) 

Po^RfR/fl + Pi{wlR - Wrl) 

^hfL/h + WLR - Wrl 
-Vl[1 - fL]/fi + WLR - Wrl 

-'^RfR/h + WLR - tlORL 
rii[l - fR]/h+WLR - WRL 



(B21) 

(B22) 

(B23) 

(B24) 
(B25) 

(B26) 

(B27) 
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A lengthy but straightforward calculation shows that 
the rate-equations results for the zero-frequency noise 
in leading and next-to-leading order are identical with 
Eqs. (B17) and (B18). 



APPENDIX C: APPROXIMATE EVALUATION 
OF COTUNNELING RATES 

We derive an approximate expression for the elastic 
cotunneling rate Wqq.^j, valid for strong electron-phonon 
coupling A ^ 1. Our starting point is the (imrcg- 
ularized) rate as obtained by Fermi's golden rule, see 
Eq. (10). This expression can be greatly simplified by 
noting that only a few terms in the q" sum give significant 
contributions, which is due to the behavior of Franck- 
Condon matrix elements for strong electron-phonon cou- 
pling. Specifically, the relevant matrix elements read 



(CI) 



exp \^—)^ + 2qhi\ + q — qhiq— 5 In 



where wc have used Sterling's approximation for the fac- 
torial in the last step. By inspection for A ^ 1, we find 
that the argument of the exponenial function is strictly 
negative and exhibits a single maximum aX q « A^. 
We obtain a reasonable approximation for the strong- 
coupling FC matrix elements by expanding the argument 
to second order in q around the maximum q = which 
leads to the Gaussian 



2-kX 



exp 



2A2 



(C2) 



Consequently, as a function of q the (squared) FC matrix 
element |Mog|^ exhibits a single peak of width ~ A cen- 
tered at g = A^. Accordingly, in Eq. (10) we only need 
to account for contributions from these summands. To 
leading order in 1/A, we can then neglect the variation of 
1 / q" within the Gaussian peak, i.e. setting q" = }? in the 
energy denominator, and noting that |e| <C (e^ — A^fio)) 
for bias and gate voltages small compared to the polaron 
shift, we obtain the approximation 



FaFfa keT + eV 



VgTb keT + eV 
TaFfc kBT + eV 



dq 



2nX 



Oq"\ 



exp 



(C3) 



{q-Xr 



2A2 



APPENDIX D: FORMALISM FOR 
CALCULATING CURRENT AND SHOT NOISE 

We review Korotkov's technique for noise calculations, 
Ref. 32, and slightly generalize his formalism in order 
to incorporate higher-order processes such as cotunnel- 
ing. Throughout this section, we will keep wording and 
formalism as general as possible to stress the wide appli- 
cability of this approach. 

We consider a system consisting of some central struc- 
ture, such as a quantum dot or molecule, and source and 
drain electrodes coupled to the central structure. Wc as- 
sume that the state space of the center is discrete and 
finite, and we denote it by 



§EE{l,2,...,iV}. 



(Dl) 



Then, the system's dynamics originates from transitions 
between pairs of these states, \i) | /), which are as- 
sociated with specific transition rates Wij > 0. Mathe- 
matically speaking, the dynamics can be classified as a 
time-dependent Markov process, and it is fully charac- 
terized by the Master equations 



d_ 
dt 



P{f, t\i) = J2 [Pik, t\i)Wkf - P{f, t\i)Wfk] (D2) 
fees 



with initial condition P{f,t = 0\i) = Sif. Here, P{f,t\i) 

denotes the conditional probability for the system to oc- 
cupy the state | / ) at time t, given that the initial state 
at time t = was \i). Note that, due to our choice of la- 
belling processes and corresponding rates with the center 
degrees of freedom only, rates Wu with identical initial 
and final states exist and need to be taken into account. 
[A relevant example is given by elastic cotunneling pro- 
cesses.] 

It is convenient to translate Eq. (D2) into a compact 
matrix notation. This can be achieved by the definitions 



Pi{t) = (P(l, t\i), P(2, t\i), . . . , P{N, t\i)f , (D3) 



and 



W 



/-Ek^iWik W21 W31 

^12 -Efe#2^2fe W32 

1^13 W23 -Efe^3^3fe ••• 



V 



(D4) 

where all /c-sums run over the state space S. Now, the 
Master equation can be summarized into 



(D5) 



with initial condition pi(t = 0) = e^, having the formal 
solution 



Pi(i)=e™*e,. 



(D6) 
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Under rather general conditions, typically satisfied in 
physical applications, one can show that for long times 
t, the probability distribution converges to a unique sta- 
tionary distribution P. Instead of performing the limit 

P = lim Pi{t) (independent of i), (D7) 

one can solve the stationary master equation 

WP = (D8) 

with the condition that the sum over all components of 
the vector P is 1 (normalization), which we write as 
trP = 1. 



1. Current 

Electron transfer between the electrodes and the cen- 
tral system gives rise to a current, whose time- average 
(J) (the stationary current) must be identical in the left 
and right junction, i.e. (7) = (II) = {Ir), due to charge 
conservation. [Here, and in the following, (•) always de- 
notes time averaging.] At this point, we may fix the sign 
of the current once and for all, defining the current to be 
positive for an electron transfer from the left to the right. 

Some caution is now in place regarding the nature of 
the transitions \i) ^ \ f) . In general, there may be sev- 
eral different processes which cause a transition between 
the central system's states | i ) and | / ) , and they differ 
from each other in the resulting leads' states. In order 
to elucidate this situation, we give one concrete exam- 
ple. Consider a system, whose central system is simply 
a single electronic level. Then, cotunneling will result 
in a transition n ^ n, i.e. the electronic occupation of 
the center remains invariant. However, there are several 
incoherent contributions to this transition: The cotun- 
neling transition can result in a transfer of one electron 
from the left to the right, or vice versa. In principle, the 
tunneling back and forth between only one lead and the 
center is possible as well. 

This example reveals that, in general, each transition 
\i) I / ) may be made up of several elementary pro- 
cesses contributing to this transition. Accordingly, the 
rates Wj/ are a sum of the rates of the corresponding 
elementary processes v, 

[Of course, the notion of elementary processes only comes 
about due to our choice of expressing all results in terms 
of the states of the central system alone.] 

Generally, not all elementary processes contribute to 
the current, and each process which does, can transfer 
a number n of electrons either in the positive or nega- 
tive direction. Taking this into account, we define the 
quantity s"^ as follows: If the elementary process v be- 
longing to the transition \ i) \ f) gives a current con- 
tribution in junction a {a = L,R), then s'^j, ^ gives the 



number of electrons transferred across junction a within 
this process, and its sign reflects whether the transfer 
occurs in the positive or negative direction. For all pro- 
cesses without current contribution, s"j^^ vanishes. Now 
the stationary current in junction a can be written as 

(4)/e = ^ = trW/P, (DIO) 

where (W^z),, =E.s,%W^^r^- 

For the analysis of the current shot noise, we are also 
interested in the explicit time-dependence of the current. 
Viewing the system dynamics as a Markov process con- 
sisting of quasi-instantaneous jumps between different 
states, we can define the times tn^'"^ {n = 1,2,..., oo) 
at which the process u belonging to the transition \i) ^ 
I / ) takes place. Then, the time-dependent current in 
junction a can formally be written as a sum of 6 func- 
tions according to 

Taking the time average of this expression and comparing 
with Eq. (DIO), we obtain the relation 

^{6it-4M)) = p,w^^;;\ (D12) 

n 

2. Noise 

Our central goal consists of the derivation of a compact 
and general expression for the current noise defined by 

/oc 
dre^"^[(4(T)4(0))-(J)2]. (D13) 
-oo 

The starting point for our considerations is the current- 
current correlator {Ia{t)Ib{t')) . When substituting 
Eq. (DIO) into the correlator, we obtain a multiple sum 
with each summand being proportional to the average of 
a product of two S functions: 

(4(t)4(0) =6^ E E E 12'k^'i'f',, (D14) 

As the set of all times tn'^''^'^ forms the time lattice of the 
Markov process, the times are considered to be pairwise 
distinct. Thus, we havetP^) ^i^f''^^ for {n,zf,i.) ^ 
(m, i' f , fJ,). Accordingly, the correlator can be separated 
into autocorrelation and cross-correlation contributions, 

{Ia{t)h{t'))a = {Ia{t)Ib{t')). + {Ia{t)W)U where 

{lamit')). = E E (D15) 

(/a(t)/&(i'))x = "non-diagonal terms". (D16) 
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We first turn to the evaluation of the autocorrelations. 
Using the fact that {S{t - - 4'^'''^)) = S{t - 

t'){S{t - t^n^'"^))), and applying relation Eq. (D12), we 
obtain 

(D17) 

The cross-correlation terms now involve only distinct 
processes at different times. These contributions can be 
dealt with in the following way. Let us assume that t>t'. 
By using I = dQ/dt, we note that 

(,.W4(0).= <f?»iM, (DI8) 

and we only get a contribution to the time average, if 
there are current-contributing processes [causing charge 
variations dQ] in jimction h during [t'.t/ + d.f/], and in 
junction a during [t,t + dt], which leads us to 

{Ia{t)h{t'))^ (D19) 

We are now in the position to give an intermediate result 
for the bracketed expression in Eq. (D13). We find 

This expression is valid for r > 0. For r < we exploit 
the symmetry 



can be expressed as 



{5UT)5hm = {5Ia{0)6h{-T)). 



(D21) 



Now, the last missing piece for an evaluation of the noise 
consists of the determination of the Fourier transform of 
the term 



Gkj{T)=e{T) \p{k,r\j)-Pk]. 



(D22) 



Noting that p{k,T\j) = {e^^)k3 and defining the matrix 
V by Vkj = Pfe, we can rewrite Eq. (D22) in matrix form 
as G(t) = 9{t) [e** — V] , and its Fourier transform is 
found to be 

/■oo 1 

Qiuj) = / dte'''^ \e^^ - V] = -(W + iojl)'^ + —V, 
^0 

(D23) 

which can directly be evaluated at all finite frequencies. 
(The w ^ will be discussed below.) Thus, the noise 



(^) 



-I- 



E Y.'kA},,PiW^pG,,{u:)Wi^^ 



(D24) 



-|- (same term with a 6)* 
= 2e^ [tr Ua6 + W6G(a;)ya -|- WaG(a;)*y6] . 
Here we have defined the vectors Ua6, yo and by 



fes f 

w. = EE4>^'^?^ 
K)^ = EE«^/..<^ 



(D25) 
(D26) 
(D27) 



In order to evaluate the zero-frequency limit a; ^ 0, we 
make use of the following properties of the relevant ma- 
trices. First, we observe that CP = 0. Second, one easily 
confirms that Vx = for any traceless vector x. Third, 
we note that any well-behaved rate-equations matrix W 
(which has a unique steady-state solution) is invertible in 
the subspace of traceless vectors. Thus, in Eq. (D24) we 
may modify the vector y^ by adding some multiple of the 
stationary distribution P, leaving the result invariant. 

We choose 



iya)j = {ya)j - {Ia)Pj/e, 



(D28) 



which makes ya a traceless vector. Consequently, the 
zero-frequency noise can be evaluated via 

Sab{<^ = 0) = 2e^ [tiuab - W6W-iy„ - w^W-iyb] . 

(D29) 



APPENDIX E: ZERO-FREQUENCY NOISE AND 
FANO FACTOR FOR THE TELEGRAPH 
REGIME 

According to the simplest situation where telegraph 
noise occurs, we take into account the transitions (1) 
I 0, ) ^ 1 1, ) and I 1, ) ^ I 0, ) due to sequential tun- 
neling, and (2) I 0, ) ^ I 0, ) and 1 1, 1 ) ^ 1 1, 1 ) due 
to cotunneling. At finite bias the transport is essentially 
unidirectional, so that we restrict to LR cotunneling, as 
well as sequential tunnehng |0,0) | 1,0) via the left 
jimction and | 1,0) — > |0,0) via the right junction. In 
order to avoid irrelevant indices, we will here denote the 
rates as W~^, W~ for sequential tunneling, and 
for cotunneling. The rate-equation matrix then reads 



W : 



'W+ w- 
w- -w+ 



(El) 
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and the resulting steady-state occupation probabilities 
for the empty and singly-occupied molecule are 



Pn = 



w- 



W+ + W- 



Pl 



w+ 



W+ + W- 



(E2) 



Note that due to the spin-degeneracy of the singly- 
occupied state, the sequential-tunneling rates obey 
W+ = 2W-, leading to Pq = 1/3 and Pi = 2/3. Ne- 
glecting sequential-tunneling rates as compared to cotun- 
neling rates where appropriate, we find that the vectors 
Uaa', Ya, and Wa are given by 



U„,„' 



w„ 



PoW° 
PiW^ 



(E3) 
(E4) 



The steady-state current is 



{I)=PoW° + PiW^ = ^W° + ^W\ (E5) 



so that Ya is obtained as 

_ w+w-{w° -W^) [ 1 

Finally, we have 



-1 



w+w-{w° -W^) ( 1 



{w+ + w-f 



-1 



For the noise, we therefore find Sll = Srh 
Srl, and the final results are 



Fps6 



w+w-{w" -W^)^ 
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This statement strictly holds for the equilibrated case 
r ^ 0. For the uncquilibratcd case r oo, the dynamics 
within avalanches involving highly excited phonon states 
may acquire small negative corrections due to cotunneling. 
More generally, the suppression applies to all sequential 
processes which leave the vibrational state unchanged or 
increase the excitation level. 



The behavior of these absorption-induced vibrational side- 
bands, which arc not specific to the case of strong clectron- 
phonon coupling, deserves further investigations, and will 
be published elsewhere. (M. C. Liiffe, J. Koch, and F. von 
Oppen, in preparation) 



